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Abstract 

The effect of solute hydrogen on the stability of vacancy clusters in hexagonal closed packed zirconium 
is investigated with an ab initio approach, including contributions of H vibrations. Atomistic simulations 
within the density functional theory evidence a strong binding of H to small vacancy clusters. The hydrogen 
effect on large vacancy loops is modeled through its interaction with the stacking faults. A thermodynamic 
modeling of H segregation on the various faults, relying on ab initio binding energies, shows that these faults 
are enriched in H, leading to a decrease of the stacking fault energies. This is consistent with the trapping 
of H by vacancy loops observed experimentally. The stronger trapping, and thus the stronger stabilization, 
is obtained for vacancy loops lying in the basal planes, i.e. the loops responsible for the breakaway growth 
observed under high irradiation dose. 
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1. Introduction 


Macroscopic properties of metals and their al¬ 
loys are known to be strongly affected by hydrogen, 
which is a common impurity found in structural 
materials. Among hydrogen effects, one can men¬ 
tion phase transformations leading to hydrides [l|, 
interactions with structural defects modifying the 
plastic activity of materials Q, the interplay with 
point defects provoking swelling |3j, or the so-called 
H embrittlement, via different possible mechanisms 


Hexagonal closed packed (hep) zirconium alloys, 
that are used as a cladding material in nuclear in¬ 
dustry, suffer an in-service hydrogen pickup due to 
the oxidation of the rods by water. Under irradi¬ 
ation, the apparent solubility limit of hydrogen in 
zirconium increases |7. [g, !). [l(| . This has been as¬ 
sociated with a trapping of hydrogen by the defects 
created by irradiation, in particular the vacancy 


clusters @. This hydrogen in solid solution influ¬ 
ences not only the mechanical properties of hep Zr 
alloys mm , but also their macroscopic structure: 
it impacts the stress free dimensional c hang e experi¬ 
enced by these alloys under irradiation [13J. Indeed, 
under irradiation, a Zr single crystal undergoes an 
elongation along the (a) axis and a shortening along 
the (c) axis, without significant volume change 14]. 
The growth strain remains small at low irradiation 
dose, whereas a breakaway growth is observed at 
higher fluence 00 - This breakaway growth at 
high doses has been correlated to the appearance of 
faulted vacancy dislocation loops lying in the basal 
planes. These vacancy clusters are called (c) loops 
because of the (c) component of their Burgers vec¬ 
tors. At low irradiation dose, where the growth 
strain remains moderate, only perfect dislocation 
loops, either of interstitial or vacancy types, with 
(a) = 1/3(1120) Burgers vector and with an habit 
plane close to the prismatic planes of the hep struc¬ 
ture, are observed. 
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Experiments have shown that an increase of the 
H content leads to an increase of irradiation growth 
in Zr alloys under neutron irradiation [l3l |. This 
is consistent with the TEM observations of Tour- 
nadre et al. 16, Et who demonstrated that the 
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amount of (c) loops formed under proton irradia¬ 
tion is notably higher when pre-hydriding the zirco¬ 
nium alloys. This suggests that hydrogen promotes 
the stability of these (c) loops compared to other 
vacancy clusters, an hypothesis also supported by 
the calorimetry and X-ray diffraction experiments 
of Vizcaino et al. These experiments showed 

that a higher annealing temperature is needed in 
zirconium alloys exposed to a high irradiation dose, 
hence when (c) vacancy loops are present, for the 
H solubility limit to reach its unirradiated value. 


Understanding the potential trapping of hydro¬ 
gen by irradiation defects as well as the hydrogen 
influence on the stability of the different vacancy 
clusters is therefore of prime importance for mod¬ 
eling then the kinetic evolution of Zr alloys under 
irradiation and the associated macroscopic behav¬ 
ior in presence of hydrogen. In this paper, we study 
the elementary interactions between H and vacancy 
clusters, going from small clusters containing only 
a few vacancies to large dislocation loops. Atom¬ 
istic simulations are a well-suited tool to address 
this question, since they give information at defect 
sizes not accessible by other techniques. In partic¬ 
ular, ab initio methods, having a high level of accu¬ 
racy and transferability, will be used all along this 
work, as they appear unavoidable to study struc¬ 
tural defects in hep Zr. A proper account of an¬ 
gular contributions is indeed necessary to correctly 
describe the stacking fault energies and the inter¬ 
actions between vacancies 18, Bdi Previous ab 
initio studies of hydrogen interaction with defects 
in zirconium mainly focused on plasticity and em¬ 
brittlement 11, 12| or on single point-defect prop¬ 


erties 21, |22[. The work presented in this paper 


extends thus these studies by modeling H interac¬ 
tion with vacancy clusters of various sizes. 


The paper is organized as follows. We first inves¬ 
tigate the interactions of hydrogen with small va¬ 
cancy clusters and with the different stacking faults 
of interest. We then model the evolution of the 
stacking fault energies with the H bulk concentra¬ 
tion. Finally, including this variation into analyt¬ 
ical laws describing the formation energies of ex¬ 
tended vacancy clusters, we discuss the influence of 
H on the relative stability of vacancy dislocation 
loops in hep Zr. 


2. Modeling method 


2.1. Ab initio simulations 

In this work, all the ab initio calculations are 
based on the Density Functional Theory (DFT), 
as implemented in the Pwscf code of the Quan¬ 
tum Espresso package [ 23 } . Calculations are per¬ 
formed in the Generalized Gradient Approximation 
with the exchange-correlation functional of Perdew, 
Burke and Ernzerhof [24]. Valence electrons are de¬ 
scribed with plane waves, using a cutoff of 28 Ry. 
The pseudo-potential approach is used to describe 
the electron-ion interaction. For Zr and H, ultrasoft 
pseudo-potentials of Vanderbilt type have been cho¬ 
sen, including 4s and 4p electrons as semicore in the 
case of Zr. The electronic density of state is broad¬ 
ened with the Methfessel-Paxton function, with a 
broadening of 0.3 eV. The integration is performed 
on a regular grid of 14 x 14 x 8 k-points for the 
primitive cell and an equivalent density of k-points 
for larger supercells. 

Atomic relaxations are performed at constant 
volume, using a conjugate gradient algorithm. For 
calculations involving point defects, we use super¬ 
cells of 5 x 5 x 4 repeated unit cells (200 atoms) 
with fully periodic boundary conditions. The elas¬ 
tic correction of Ref. [25| is applied, so as to remove 
the spurious interaction energy between the defect 
and its periodic images. |Appendix A| demonstrates 
that the obtained corrected energies are well con¬ 
verged with respect to the supercell size. 


2.2. Vibrations 

To consider the contribution of vibrations in 
the computed energies, we use the harmonic ap¬ 
proximation [ 26 }, which validity is assessed in 
| Appendix B| for both a single H atom and an H 
interacting with a vacancy in an hep Zr matrix. 
Within this approximation, the vibrational free en¬ 
ergy is 


F vib (T) = kT^2 In 

i 

with the sum running over all vibration modes of 
pulsation u>i . We further make the assumption that 
the vibration modes of the Zr atoms are not affected 
by the presence of hydrogen and consider only the 
three vibration modes of the H atom. This is jus¬ 
tified by the small mass and the small radius of 
hydrogen compared to zirconium. 
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To obtain the H vibration frequencies for a given 
atomic configuration, we rigidly displace the H 
atom from its equilibrium position in six differ¬ 
ent directions, using four intermediate positions 
for each direction with a maximal displacement of 
0.04 A. We then make a least-square fitting of the 
calculated energy variations to obtain the 3x3 Hes¬ 
sian matrix, from which the eigenfrequencies are 
obtained by diagonalization. 


2.3. Validation of the modeling approach 

Our choice of cutoffs, k-mesh, GGA functional 
for the exchange correlation and pseudo-potential 
for zirconium have already been validated both on 
the hep bulk and on vacancy cluster properties in 
previous studies 0,0 Si For H, the validation 


of the chosen pseudo-potential consists first in com¬ 
paring with experimental data the intrinsic prop¬ 
erties of the H 2 molecule, such as the equilibrium 
distance du 2 , the vibrational frequency zzh 2 and the 
dissociation energy E^ ss (see Table [TJ. The dis¬ 
sociation energy includes both the binding and vi¬ 
brational \hv h 2 contributions. The vibrational fre¬ 
quency vn 2 is calculated according to the harmonic 
oscillator approximation: 



where k is the bond stiffness obtained by ab initio 
calculations and ji = tor/2 is the reduced mass of 
the H 2 molecule. The obtained quantities are all in 
good agreement with the experimental ones, within 
less than 5%. This validates our choice of pseudo¬ 
potential for the calculations involving hydrogen. 

As we are interested in the interplay between hy¬ 
drogen solute and vacancies in an hep zirconium 
matrix, we also check the insertion of the H atom in 
substitutional position and in different interstitial 
sites: the tetrahedral (T) and the octahedral (O) 
sites. Our DFT calculations show that the substi¬ 
tutional position for the H atom, corresponding ac¬ 
tually to an H atom inside a vacancy, is metastable 
and has a higher energy than the interstitial posi¬ 
tions. More details are given in the following section 
of this paper. For the interstitial sites, we calculate 
the solution energy, defined as follows: 

-Er 1 = -EzrH — -Ezr — ^ (3) 

where Ez r refers to the energy of the Zr bulk su¬ 
percell, F/r 2 to the energy of a H 2 molecule (in¬ 
cluding the vibrational contribution) and F?zi-h to 


the energy of the supercell containing an H atom 
located in the different possible sites. The calcu¬ 
lated solution energy, without the contribution of 
H vibration in the Zr matrix, is lower for H in the 
T site than for H in the O site, with an energy 
difference A E°/ T = 60meV, in good agreement 
with other ab initio calculations although variations 
are observed between different studies (72meV in 
0, 80meV in 0, 86meV in 0, 40meV in 0, 
and GlmeY in [22[). Neutron diffraction experi- 


31] have shown that the H atoms occupy 
The correct interstitial site 


ments 
the T sites in hep Zr. 


is thus predicted by ab initio calculations and our 
value of solution energy Flff 1 is very close to the 
experimental one [0 (see Table [T|. But this is 
not true anymore when including the H vibration 
energy in the calculation. With this contribution, 
the energy difference between the T and O inser¬ 
tion sites becomes A E°/ T = — 26meV at OK and 
A E°/ t = — 47.meV at 600 K, thus corresponding 
to the O site being more stable than the T site, 
with an increasing stability with the temperature. 
As pointed by Christensen et al. 0] and shown in 
| Appendix B.l[ the energy difference between these 
two interstitial sites is too small to be able to con¬ 
clude on the preferential insertion site from ab ini¬ 
tio calculations without a much more elaborated 
treatment of vibrations. Consequently, both T and 
O insertions sites are considered in the following. 
Although it may not be the most stable when vi¬ 
brations of H are considered, we take the T insertion 
site as the H reference state in the hep Zr matrix 
for the various calculations of binding energies. 

We also verify that the insertion of H 2 in Zr is 
not favorable: it is unstable inside a T site, and 
metastable inside an O site, but with a negative 
binding energy corresponding to a strong repulsion 
(F1 r 2 = —0.70 eV). This shows that hydrogen at 
the molecular state does not need to be considered 
at low H concentrations. 


3. Hydrogen and small vacancy clusters 

We first examine the interaction between H and 
small vacancy clusters, denoted as V ra , with n the 
number of vacancies. We define the binding energy 
between H and the vacancy cluster as: 

Ar Vii = -EzrVn + -®ZrH ~ -EzrHV„ — -E Zr (4) 

where i?zrV„, -EzrH, FlzrHVn and -®Zr are the ener¬ 
gies of the same simulation cell containing respec¬ 
tively the vacancy cluster, an H atom in T site, 
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Table 1: Properties of the H 2 molecule (equilibrium distance djj 2 , vibration frequency //h 2 and dissociation energy E^ s ) 
and solution energies E ^ ol of H in the hep Zr matrix. Values in parenthesis include zero point vibration energy of H in the Zr 
matrix. 



This work 

A b initio 
Udagawa [12] 

i 

Domain [28] 

Lumley [29] 

Exp. [27] 

dn 2 (A) 

0.75 

0.75 

0.76 

0.75 

0.74 

^ (THz) 

130.1 

- 

134.4 

129.8 

130.8 

(eV) 

Elf (eV) T site 

4.25 

4.50 

4.27 

4.53 

4.48 

-0.609 (-0.387) 

-0.52 

-0.604 

-0.60 

-0.66 

Eff (eV) O site 

-0.549 (-0-414) 

-0.44 

-0.532 

-0.56 

- 


the complex HV„ and no defect. A positive value 
of -Ehv indicates an attractive interaction between 
the H atom and the vacancy cluster. For a given 
cluster V n , different possibilities exist to insert the 
H atom. It can be located either in a substitu¬ 
tional or in an interstitial site (T and O, here), and 
many of these sites can be found close to the va¬ 
cancy cluster. All of them should be considered, 
but the amount of different configurations to ex¬ 
plore strongly increases with the cluster size, mak¬ 
ing an exhaustive study out of reach. An approach 
based on the understanding of the elementary inter¬ 
actions between H and vacancies therefore appears 
necessary, so as to reduce the configuration space 
and to focus on important clusters. 


[ 7 ] from a kinetic modeling of his recovery exper¬ 
iments performed after deuterium implantation in 
irradiated zirconium. 
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3.1. H-V interaction 

We start by studying the interaction between H 
and a single vacancy. We first check the direct in¬ 
sertion of H into a vacancy. It is metastable, but 
with a strongly repulsive interaction: we obtain 
v = —1.16 eV for this position, a value close 


to the results of Domain et al. [28|. So we do not 


consider anymore this substitutional position in the 
following. 

We then examine both types of interstitial sites, 
T and O, for the insertion of hydrogen in the vicin¬ 
ity of the vacancy. Different configurations are in¬ 
vestigated, where the hydrogen atom is placed in 
the successive neighboring shells of the vacancy. 
The configurations corresponding to the different 
nearest neighbor positions are displayed in Fig. [I] 
The resulting binding energies between H and V are 
provided as a function of their separation distance 
in Fig. [2] The maximal binding energy is 0.21 eV 
without H vibration and 0.24 eV with the inclusion 
of zero point energy, in good agreement with the ex¬ 
perimental value 0.28 ± 0.06 eV obtained by Lewis 


Q.O 

/•'■■■. /o\ 

O Q O 

□ , / 

O O 

Configuration c 


Q.O 

■ O \ 

P O Q 
. O ■■!■■■ O ■■ 

o □ o 

Configuration a’ 


Figure 1: Projection in the basal plane of the different H- 
V pairs investigated. The successive neighbor shells of an H 
atom lying either in a T or O interstitial site are respectively 
denoted a to g and a’ to d’. The white spheres represent the 
Zr atoms at z = 0, the grey spheres the Zr atoms at £ = c/2 
and the squares symbolize the vacancy. The purple spheres 
represent the H atom in an O site at z — c/4, the green 
spheres the H atom in a T site at 2 : = c/8 (at z = —c/8 for 
the small symbol used in configuration c). 

We first analyse the case where H is in a T site. 
The two non-equivalent first neighbor positions lead 
to different binding energies: the value for the con¬ 
figuration b is twice the one of the configuration 
a. Accounting for the directional aspect of the H-V 
interaction therefore appears to be important. The 
c/a ratio of the hep Zr structure is 1.601 in our 
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distance (A) 

(a) without H vibration 



2 4 6 


distance (A) 

(b) with H vibration 

Figure 2: Binding energy between H and a vacancy versus 
their separation distance, (a) without and (b) with H vibra¬ 
tion free energy. Both O and T insertion sites are considered 
for H. Label sites a to g and a’ to d’ are detailed in Fig. ^ 

ab initio approach 0 = this is a lower value than 
the ideal one (-^/8/3). Imposing an ideal c/a ratio, 
we calculate again these two binding energies. We 
obtain in that case E\ v = 0.04 eV for the configu¬ 
ration a and E\ v = 0.20 eV for the configuration 
b. The two interactions remain different even when 
applying an ideal c/a ratio. The in-plane interac¬ 
tion (configuration b) is rather unaffected by the 
c/a ratio, while the out-of-plane interaction of con¬ 
figuration a is very sensitive to its value. It goes 
from attractive to zero when increasing c/a: con¬ 
sequently, we expect this binding energy between 
H and V to be strongly affected when applying a 
strain. Configuration c shows an attractive binding 
of slightly larger magnitude than the one of config¬ 
uration a. For larger distances (configurations d to 


g), the binding energy rapidly goes to zero, show¬ 
ing that long range interactions do not exist in this 
case. 

We finally look at the insertion of H into the O 
site. Its binding energy with the vacancy is attrac¬ 
tive only for the nearest neighbor configuration, de¬ 
noted a’ here. Its magnitude is close to the mag¬ 
nitude of the most attractive configuration for H 
in T site (configuration b). This is surprising, as 
in pure zirconium, H is located in the tetrahedral 
sites: the presence of the vacancy hence increases 
the stability of the octahedral site. This is an addi¬ 
tional reason to always consider this insertion site 
in the rest of our study. For the configurations b’ to 
d’, the binding energy decreases to reach the energy 
difference between O and T sites in bulk Zr, thus 
indicating that H does not interact anymore with 
the vacancy. 

The main change caused by the vibrational free 
energy is a constant shift of the binding energies for 
the configurations with an H atom in an O site (Fig. 
m This merely corresponds to the increase of sta¬ 
bility observed for the O site in bulk Zr. Apart from 
this constant shift, H vibrations slightly modify the 
interaction between the H and the vacancy when 
they are first nearest neighbours: it increases the 
binding when the H atom is in a T site (configura¬ 
tions a and b) and decreases it for the O site (con¬ 
figuration a’). When the H atom and the vacancy 
are separated by a larger distance, the H vibrations 
does not change their binding. 

This analysis of the H-V interaction suggests that 
only the very short ranged interactions need to be 
considered for the stability of the complexes be¬ 
tween hydrogen and vacancies. This point is dis¬ 
cussed in the following. 

3.2. Position of H close to a vacancy cluster 

Vacancies are clustering under irradiation in zir¬ 
conium alloys, and we have shown that H and V 
interact at very short distances. In order to under¬ 
stand where the H atom is located with respect to 
a preexisting vacancy cluster, we consider now the 
accumulation of vacancies around a single H atom. 
To this aim, we calculate the binding energies, as 
defined in Eq. QJ of H with clusters of n vacancies, 
where n vacancies surround the H atom in near¬ 
est neighbor positions, i.e. vacancy positions corre¬ 
sponding to configurations a, b and a’ for the H-V 
pair (Fig. [T]). Again, both T and O interstitial sites 
are considered. For each number n of vacancies, 
the different possible configurations are tested (see 
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Tab. [2j, and the most stable one is retained. The 
contribution of H vibrations to the binding energy 
is considered only for these most stable clusters. 
The resulting binding energies versus the number 
of vacancies surrounding the H atom are displayed 
in Fig. [3] H vibrations generally increase the sta¬ 
bility of the HV„ clusters. 

Table 2: Different possible configurations of the H-V TO com¬ 
plexes involving n nearest neighbor bonds between H and 
V. The H atom is located either in T or O site. Values in 
parenthesis include zero point vibration energy of H. 
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For H in T site, the binding energy decreases with 
the number n of nearest neighbor vacancies, except 
at high enough temperature when the H vibrations 
are considered. The nearest neighbor interactions 
are non additive and the highest binding energy is 
obtained for the HV complex. In addition, some 
configurations, involving 2, 3 and 4 nearest neigh¬ 
bor vacancies, are unstable. This suggests a more 
favorable insertion of the H atom at the boundary 
of a vacancy cluster, rather than inside the cluster. 

The situation is quite different when H is inserted 
into the O site. Up to 3 vacancies, the binding 
energy is partially additive: adding a new nearest 
neighbor vacancy to the HV n cluster stabilizes it a 
bit more. For n > 3, the binding energy rapidly 
decreases, becoming repulsive for n = 5. For n = 6, 



number of nearest neighbor vacancies 
(a) without H vibration 



number of nearest neighbor vacancies 
(b) with H vibration 

Figure 3: Binding energies between H and an increasing 
number of vacancies in nearest neighbor position, (a) with¬ 
out and (b) with H vibration free energy/ s f , 5c H atom is 
located either in a T or O site. 
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the cluster corresponds to an H atom inside a small 
cavity. As a consequence, the HV 6 cluster is unsta¬ 
ble, with negative frequencies for the H vibration. 
Again, this suggests that the insertion inside a va¬ 
cancy cluster is not favorable. Interestingly, the 
configurations where H occupies an O site with 2 
or 3 first nearest neighbor vacancies lead to higher 
binding energies than any configuration where H is 
in a T site. This was not expected since H is lo¬ 
cated in the T site in pure bulk hep zirconium. A 
similar behavior for H close to vacancy clusters has 
already been observed recently in bcc iron [32| . 

3.3. Interaction with small vacancy clusters 

In light of the two preceding sections, a very sim¬ 
ple approach can be proposed so as to reduce the 
configuration space to explore for finding the most 
stable insertion sites of H close to a vacancy cluster. 
We have identified three preferential insertion sites 
for H : T site with one vacancy in nearest neighbor 
position, and O site with 2 or 3 nearest neighbor 
vacancies. Only those sites will be considered in 
the rest of our study. The magnitude of the H- 
V interaction rapidly decreases with the distance: 
we therefore suppose that only the nearest neigh¬ 
bor interactions are important, and that the longer 
ranged terms (g?hv > 3.5 A) do not influence the 
stability of H close to a vacancy cluster. Conse¬ 
quently, we assume that the binding energy between 
H and the vacancy cluster mainly depends on the 
nature and the number of nearest neighbor interac¬ 
tions between H and the vacancies. Only the config¬ 
urations that maximize the binding energy accord¬ 
ing to this assumption are considered and relaxed 
with ab initio calculations. We first check the va¬ 
lidity of this assumption, and we investigate then 
the interaction of H with vacancy clusters of differ¬ 
ent types. As H vibrations do not strongly modify 
the H interaction with vacancy clusters, we do not 
include them in our calculation of binding energies 
in this section. 

3.3.1. Validation of the approach 

We test the validity of this approach on two dif¬ 
ferent clusters containing four vacancies: one sit¬ 
ting in the basal plane and the other one sitting 
in the prismatic (1010) plane. These two clusters 
are the most stable plane clusters of 4 vacancies 
in hep Zr [§o|. We compute the binding energy of 
these clusters with H for different positions of the H 
atom involving p first nearest neighbor interactions 


(1 < p < 3) and compare with the binding energy 
of H with a cluster containing only p vacancies, and 
where the same p interactions are involved. T and 
O sites are explored, and the calculated binding en¬ 
ergies are displayed in Table [3] Calculations for this 
validation test are performed in a smaller 4x4x3 
supercell containing 96 lattice sites. 

Table 3: Comparison of binding energy of H with V 4 clusters 
where p nearest neighbor H-V interactions are present, with 
the binding energy of H with V p clusters where the same p 
nearest neighbor H-V interactions are present. Both T and 
O insertion sites are explored. 
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When H is inserted into a T site, assuming that 
the number and the nature of nearest neighbor va¬ 
cancies surrounding the H atom determines the sta¬ 
bility of the cluster provides a rather good approxi- 
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mation. First, the instabilities are well predicted by 
this simple approach (Tab. [3]), and the final configu¬ 
rations have the same structure after atomic relax¬ 
ations. Second, the differences between and 

the corresponding f?HV p are smaller than 0.07eV. 
For H in O site, the agreement between ^hv 4 an d 
is even better (Tab. [3J. The obtained val¬ 
ues for the clusters HV 4 and HV P with p near¬ 
est neighbor interactions are very close, within less 
than 0.02 eV in this case. 

This study confirms that vacancies farther than 
nearest neighbor positions have only a weak influ¬ 
ence on the binding energy value. This latter is 
mainly fixed by the number and the position of va¬ 
cancies that are nearest neighbors of the H atom. 
Again, we note on the HV 4 clusters that the bind¬ 
ing energies for H in O site with p = 2 or 3 are 
higher than the binding energies for H in T site. 


3.3.2. Influence of H on small clusters 
We finally examine the stability properties of the 
HV„ clusters for 2 < n < 7. We separate vacancy 
clusters into three different groups : 


• basal clusters, where all vacancies are lying in 
the same basal plane. These clusters can be 
seen as precursors of (c) loops. 


• prismatic clusters, where all vacancies are lying 
in the same prismatic corrugated plane. These 
clusters can be seen as precursors of (a) loops. 

• 3D clusters (precursors of cavities). 


In the absence of hydrogen, we have shown in 
a previous study [20] that the 3D small vacancy 
clusters are more stable than the plane vacancy 
clusters, until 7 vacancies. As stated in our intro¬ 
duction, under irradiation, the amount of (c) va¬ 
cancy loops increases when zirconium samples are 
pre-hydrided 16, Ezi- This could suggest that H 


in solid solution affects the relative stability of the 
different types of vacancy clusters. Studying the 
binding properties of H with small V„ clusters of 
different types should give some insight about the 
validity of such an assumption. 

We proceed as follows. The chosen vacancy clus¬ 
ters V„, with n = 2 to 7, correspond to the most 
stable ones, previously found in pure Zr [2Cj. For 
each cluster, the previously identified insertion sites 
for H are tested. The most stable one is retained 
and we plot the binding energies versus the clus¬ 
ter size for each group of cluster (see Fig. [4]). The 



Figure 4: Maximum binding energies of H with vacancy clus¬ 
ters containing between 2 and 7 vacancies. Three types of 
clusters are considered: basal, prismatic and 3D clusters. 


resulting binding energies are all positive, and are 
included in a range from 0.35 to 0.48 eV. 

No clear preference is seen for any type of va¬ 
cancy cluster. This is not surprising, since, as es¬ 
tablished above, H mainly interacts with its first 
nearest neighbor vacancies. H therefore does not 
have any discriminating effect on the small vacancy 
clusters. 


4. Hydrogen interaction with stacking faults 

In hep Zr, vacancy clustering under irradiation 
leads to planar extended defects: (a) dislocation 
loops, lying in planes close to the prismatic planes, 
and (c) component loops, lying in the basal planes. 
Cavities are hardly ever observed, in agreement 
with their lower stability [ 20 | . and will not be con¬ 
sidered in the following. The (a) loops are the most 
stable defects {2(|. They are faulted at small sizes 
and perfect at large sizes. The (c) loops on the 
other hand are always faulted, with an extrinsic 
basal fault at small sizes and an intrinsic Ii fault 
at large size. These (c) loops are believed to play a 
crucial role in the breakaway growth at high irradi¬ 
ation dose. 

In order to determine the influence of H solutes 
on vacancy dislocation loops, we perform a careful 
study of their binding energies with the different 
stacking faults involved in the extended defects of 
interest. We assume that this interaction with the 
stacking fault is the main H contribution on the 
energetics of vacancy loops. 





4-1- Prismatic stacking fault 


When removing a vacancy platelet in a corru¬ 
gated {1010} plane, a prismatic stacking fault is 
formed [§o[. This stacking fault is associated with 
the formation of faulted (a) loops of Burgers vec¬ 
tor b = 1/2 (1010). The unfaulting of the vacancy 
loop occurs by a 1/6(1210) shearing of the fault 
plane, leading to a perfect (a) loop of Burgers vector 
b = 1/3 (2110). We investigate here the interaction 
of an H atom with this intrinsic stacking fault. 

Calculations are performed in a simulation box 
duplicated 4 times in the [1010] direction, 2 times 
in the [0001] direction, and 4 times in the [1210] 
direction, i.e. in the direction perpendicular to the 
fault plane. This corresponds to a total of 128 Zr 
atoms. The use of this simulation box to model a 
prismatic stacking fault has been already validated 
[191 ]. It is also large enough to model an isolated 
H atom: with such a simulation box, the energy 
difference between the O and T insertion sites is 
AE°/ t = 57 meV in the perfect crystal, in good 
agreement with the value obtained in conventional 
supercells (Tab. [TJ . 

Different insertion positions for the H atom, ly¬ 
ing either in a T or an O interstitial site at vary¬ 
ing distances from the stacking fault, are considered 
(Fig. [5}i). Ab initio calculations lead to a binding 
of the hydrogen with the prismatic stacking fault 
(Fig. [5jo and c), with a maximal binding energy 
of 0.1 eV. The most attractive interaction is ob¬ 
tained when the H atom is one plane apart from the 
stacking fault into a T site. This attractive inter¬ 
action is in agreement with previous ab initio stud¬ 
ies, although our maximal binding energy is slightly 
lower than the one obtained by Domain et al. [11] 
(0.14eV) or by Udagawa et al. 12| (0.16eV). The 
binding energy rapidly decays when the distance 
between H and the stacking fault increases: it goes 
to zero for the T site and to —A E°/ T for the O 
site. The only octahedral site leading to an attrac¬ 
tive interaction with the stacking fault is the closest 
to the fault plane. This octahedral site did not ex¬ 
ist before the introduction of the stacking fault and 
is created by the shearing of the lattice j33|. The 
position corresponding exactly to the octahedron 
position in this site is unstable. For each such oc¬ 
tahedral site, two stable positions inside each half- 
octahedron exist, where the H atom is linked to five 
Zr first nearest neighbors (Fig. (5ji). For all other 
octahedral sites, as well as all the tetrahedral sites, 
the relaxed position of the H atom is close to the 



[ 1210 ] 

(a) H insertion sites 



distance (A) 

(b) Binding energy without H vibration 



0 2 4 6 8 10 

distance (A) 

(c) Binding energy with H vibration 


Figure 5: Interaction between H and a prismatic stacking 
fault, (a) Projected positions in the [1210] plane of the in¬ 
vestigated T and O sites for the insertion of H in the vicin¬ 
ity of the stacking fault (vertical dashed line). Zr atoms are 
shown by white and grey circles, depending on their position 
along the [1210] direction in the faulted crystal, (b) and (c) 
Binding energies ^H_'y between H and the prismatic stack¬ 
ing fault, versus their initial separation distance, calculated 
(b) without and (c) with H vibration free energy. 
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ideal position. 

The same energy shift as in the bulk, correspond¬ 
ing to a stabilization by vibrations of the O sites 
compared to the T sites, is observed in presence of 
the stacking fault (Fig. [5jc) . Apart from this energy 
shift, the H vibrations decrease the H interaction 
with the stacking fault for the O site, whereas it 
increases this interaction for the T sites. 


4-2. Basal stacking faults 

Condensation of vacancies in a basal plane re¬ 
sults in the creation of a dislocation loop of Burg¬ 
ers vector b\ = 1/2 [0001]. This corresponds to the 
removal of a platelet of one atomic layer in the per¬ 
fect stacking BAB ABA of basal planes and leads 
to the formation of a highly energetic stacking se¬ 
quence, BAB.BAB A. The stacking then evolves so 
as to lower the energy of the vacancy loop by cre¬ 
ating one of the two following stacking faults [34| : 
an extrinsic fault E, with a B ABC ABA stacking 
sequence preserving the same Burgers vector, or an 
intrinsic fault Ii, with a BABCBCB stacking se¬ 
quence leading to a dislocation loop with Burgers 
vector 62 = 1/6 (2023). We investigate here the in¬ 
teraction of H with these two basal stacking faults, 
which gives insights on the influence of this impu¬ 
rity on the stability of (c) loops. 

Ab initio calculations are performed on simula¬ 
tion boxes duplicated 3 times in the [2110] direc¬ 
tion and 3 times in the [1210] direction. To build 
the basal stacking faults, one needs to remove one 
atomic layer normal to the [ 0001 ] direction, thus 
leading to a total of 2 n z — 1 atomic layers, if n z 
is the number of duplicated unit cells in the [ 0001 ] 
direction. Our simulations are carried out using 
n z = 6 . This box size, which corresponds to a total 
of 99 Zr atoms, has already been validated for both 
basal stacking faults in a previous study [ 20 ], and 
is also large enough to simulate an isolated H atom. 
The binding energy of H to the basal stacking fault 
then reads: 


Ei 


H —7 


rp2n z — l rp2n z —1 

^ZrH I" -^Zr+7 


rp2n z — 1 ta 271 z —1 

^ZrH+7 -^Zr 


(5) 


£z"h _1 is not accessible directly, as a periodic simu¬ 
lation box with 2 n z — 1 atomic layers automatically 
includes a stacking fault. We thus make the follow¬ 
ing approximation: 


rp2n z — l _ rp2n z — 1 ^ -rp2n z _ rp2n z 

-^ZrH ^Zr — -^ZrH ^Zr 5 


( 6 ) 


which allows us to obtain the binding energies of H 
with the basal stacking faults. 

As for the prismatic fault, our ab initio calcula¬ 
tions lead to a binding of H with basal faults, but 
with a higher maximal value of 0.16 eV for the two 
basal faults. The strongest bindings are obtained 
when H is close to the stacking fault. This is true 
for both T and O sites, and for both Ii and E faults. 
For the Ii fault, the T sites are always more favor¬ 
able than the O sites to insert the H atom. A higher 
number of attractive positions in the vicinity of the 
stacking fault are found for the insertion in T sites. 
The binding energies rapidly decrease with the dis¬ 
tance, and become null or equal to —A E 0//T for 
distances larger than 6 A. For the E fault, the most 
attractive interaction is found for the O site. A 
largest number of attractive positions can be found 
in T sites than in O sites. Finally, interactions of 
H are longer ranged with the E fault than with the 
Ii one, as non-interacting behaviors are recovered 
only for separation distances larger than 9 A. This 
might be related to the fact that in the E stacking 
sequence, more atomic layers have an environment 
different from the hep one than in the Ii stacking 
sequence. The effect of the H vibration is the same 
as for the prismatic stacking fault, i.e. mainly a 
shift of the binding energies for the O sites. 


5. Hydrogen and vacancy loops 

To simulate the properties of extended defects 
such as dislocation loops, one cannot rely directly 
on ab initio methods, because of the limited num¬ 
ber of atoms, only a few hundreds, which can be 
modeled with such simulations. In a previous study 
[io| . we adopted a continuous description of large 
vacancy clusters. For dislocation loops, the clus¬ 
ter energies were obtained through a line tension 
model, incorporating the stacking fault energy for 
faulted loops. After validating the approach with 
large scale atomistic simulations relying on an em¬ 
pirical potential, the main parameters of the model 
were calibrated on ab initio data. 

The influence of hydrogen on the stability of dis¬ 
location loops can be introduced in a very simple 
way within this approach, by considering the mod¬ 
ification of the stacking fault energies arising from 
hydrogen segregation [35[. Here, we first model the 
segregation of hydrogen on the different stacking 
faults, as well as the evolution of the stacking fault 
energies with the H content. We then discuss the 
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(a) H insertion sites (Ii fault) 




0 2 4 6 8 10 

distance (A) 

(c) Binding energy without H vibration (Ii fault) 



distance (A) 

(d) Binding energy without H vibration (E fault) 



0 2 4 6 8 10 

distance (A) 

(e) Binding energy with H vibration (Ii fault) 



0 2 4 6 8 10 

distance (A) 

(f) Binding energy with H vibration (E fault) 


Figure 6: Interaction between an H atom and basal stacking faults. The different T and O interstitial sites, where an H atom 
has been inserted, are sketched in (a) for the Ii and in (b) for the E basal stacking fault. The corresponding binding energies 
E^_ are respectively shown in (c) and (d) without H vibration free energy, and in (e) and (f) with H vibration free energy. 


effect of H on the stability of the different vacancy 
loops formed under irradiation in hep Zr. 


5.1. Segregation profiles 

Following common approaches for solute segrega¬ 
tion at a surface or an interface [Ullzl , we make 
use of the previously calculated binding energies 
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Eft to model the equilibrium concentration pro¬ 
files of H in the vicinity of each stacking fault. 

We consider a bicrystal consisted of two semi in¬ 
finite bulk spaces separated by an interface corre¬ 
sponding to the stacking fault. The different inser¬ 
tion sites for H, either T or O sites, are sitting in 
planes p located at different distances d p from the 
stacking fault. We assume that the system is large 
enough so that each side of the stacking fault acts 
as a reservoir for hydrogen. The hydrogen chemi¬ 
cal potential pn is then fixed at equilibrium. We 
assume that the concentration c p G [0,1] of each 
plane is homogeneous at equilibrium and that the 
configurational entropy is equal to the one of an 
ideal solid solution. We neglect the interaction be¬ 
tween H atoms and only the H interaction with the 
stacking fault is considered, a valid assumption for 
not too concentrated solid solutions. Minimizing 
the grand potential of the system with respect to 
the average H occupation c p of each plane p, the 
equilibrium distribution obey the equations 

-^h 1 — Eft_ 7 (d p ) + kT In —^— = m , (7) 

where E g? 1 is the hydrogen solution energy in the 
perfect zirconium matrix for the T insertion sites. 
These equations are valid for any plane of the 
bicrystal, as each plane, defined by its distance d p 
to the fault plane, has the same density of inser¬ 
tion sites. This is true both for the prismatic (Fig. 
[5^) and the basal faults (Fig. [HJi and b), the inser¬ 
tion sites being either of the T or of the O type. 
Far from the fault plane, concentrations and ener¬ 
gies converge to their bulk values: c p = and 
£& ((*«,) = 0 for planes corresponding to T inser¬ 

tion sites, c p = Cq and Eft (doo) = —A E°/ T for 
O planes. This allows us to eliminate /in in Eq. |7] 
and to obtain for each plane p: 


in the low concentration regime (xh <C 1). We 
finally obtain the concentration profile as a function 
of the nominal concentration: 


xn exp 


E&-^d p ) 

kT 


2 + exp (- A fc^ /T ) + x H exp ( 


(d p ) 


kT 


r 

(ii) 



(a) without H vibration 



distance (A) 
(b) with H vibration 


Cp 

1 Cp 



exp 



( 8 ) 


There are two interstitial site of T type and one of 
O type per Zr atom. Their bulk concentrations is 
therefore linked to the nominal hydrogen concen¬ 
tration xr trough the equations: 


= 


Co = 


2 + exp (-^) IH ' 

(-^) 


exp i 


! + exp(-^r 


(9) 

( 10 ) 


Figure 7: Concentration profiles of hydrogen at T = 600 K 
and for a nominal concentration ith = 0.01, plotted as a 
function of the separation distance from the prismatic, basal 
Ii and basal E stacking faults, (a) neglecting or (b) consid¬ 
ering H vibration free energy. Small triangles refer to the T 
sites and small diamonds to the O sites. The continuous and 
dashed horizontal lines correspond to the concentrations of 
the, respectively, T and O sites in the bulk. 

Fig. [T] displays the concentrations c p for the pris¬ 
matic, basal Ii and basal E faults as a function 
of the separation distance from the stacking fault. 
The segregation profiles are given for T = 600 K, 
the service temperature of the pressurized water re¬ 
actors, and for a bulk nominal concentration xh = 
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0.01. This corresponds to an upper limit of the 
H content that can be reached, as the solubility 
limit, at this temperature, is 0.008 in pure Zr [381 ] 
and 0.009 in Zircaloy-4 [9j, a zirconium alloy cur¬ 
rently used in the nuclear industry. All faults lead 
to a strong increase of the H concentration in their 
vicinity. As these faults correspond to the ones met 
in vacancy loops, this agrees with experiments sug¬ 
gesting a trapping of H by these loops 000 . We 
also notice that the basal faults concentrate more H 
atoms than the prismatic fault, in particular the E 
fault which can capture hydrogen up to the fourth 
layer from the stacking faults. The H trapping by 
basal faults is therefore stronger, explaining why 
annealing at higher temperatures is needed to re¬ 
cover the unirradiated solubility limits in Zr sam¬ 
ples where vacancy basal loops have been formed 
after an irradiation at a high dose 0- 

H vibrations do not change these general trends 
between the different stacking faults, although they 
modify the segregation profiles. The inclusion of 
the vibration free energy in the binding energies 
leads to an increase of the enrichment for the O 
sites and a decrease for the T sites (Fig. [7]). 


5.2. Variation of stacking fault energies 


As first pointed out by Suzuki Q, solute seg¬ 
regation modifies the stacking fault energies. We 
thus need to deduce from our segregation model, 
the evolution of the basal and prismatic fault ener¬ 
gies in the presence of hydrogen. As the system is 
open with respect to H, we define the stacking fault 
energy 7 ({c p }) directly from the grand potential: 


7 ({ c p}) = (12) 

*0 

where So is the area of the dividing interface for 
one insertion site per plane: So = a 2 V 3/2 for the 
basal faults and So = ac/2 for the prismatic fault. 
flf and f \ are the grand potentials of respectively 
the faulted crystal with the segregation profile and 
the perfect crystal with the homogeneous solid solu¬ 
tion. The thermodynamic model used in the previ¬ 
ous section to obtain the segregation profiles leads 


to the following expressions: 

^4h) =y; {4 [Eg 1 - Mh] 

p,t 

+kT [4 In (4) + (1 — 4) l n (1 — 4)] } 

+ E(4 [Eg + AE°/ T - m 
p, o 

+ kT [4 In (cq) + (1 — Cq) In (1 — 4)] }, 

(13) 

%({c P }) =S 0 7o + E c p W - E u-y{d P ) ~ m] 

P 

+ kT i c p ln °p + (i — c p ) ln ( l _ °p )] i 

p 

(14) 


with 70 the stacking fault energy without hydrogen. 
The summation over the planes in Eq. [13] has been 
partitioned between the planes containing T and 
O interstitial sites. Inserting Eqs. [13] and [TJ] into 
Eq. [T7] and eliminating the chemical potential /rn 
using Eq. 0 we obtain the stacking fault energy in 
presence of the segregation profile: 


7({c P }) =7o + 


kT- 


p, t 


ln 


1 Cp 

i — 4 


+ 


kT 

5b 


5 > 



(15) 


The same result is obtained using the macroscopic 
thermodynamic approach summarized by Hirth 
0. 


Using Eq. [15] we show in Fig. [ 8 ] the evolution 
at T = 600 K of the two basal and the prismatic 
stacking faults, when the nominal hydrogen content 
varies in a concentration range Xn < 0.01. All the 
stacking fault energies decrease when the hydrogen 
content increases. The decay is small for the prism 
fault, slightly higher for the U fault and fast for the 
E fault. The consideration of H vibration only af¬ 
fects the E fault: it leads to a stronger decrease of 
the fault energy with the H content. In the range 
of nominal concentrations allowed for hydrogen at 
this temperature, the relative order of the differ¬ 
ent fault energies is not modified, except close to 
the solubility limit where the basal E and the pris¬ 
matic stacking faults have almost the same energy, 
with the basal E becoming slightly most stable with 
the inclusion of vibration free energies in the ther¬ 
modynamic modeling. 
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Figure 8: Evolution of the basal and prismatic stacking fault 
energies with the bulk concentration :ch of hydrogen calcu¬ 
lated at 600 K. The thick dashed lines and the thin con¬ 
tinuous lines correspond to the data respectively with and 
without H vibration free energy. The continuous and dashed 
vertical lines correspond to H solubility at this temperature 
respectively in pure Zr and in Zircaloy-4. 


5.3. Stability of vacancy loops 

We now study how hydrogen may modify the sta¬ 
bility of the different vacancy loops, using the con¬ 
tinuous model we previously developed to describe 
the energetics of vacancy clusters in pure zirconium 
[20I ]. Within this model, the formation energy of a 
loop containing n vacancies is given by: 

E[ oop (n) = 7uRi 2 7n + 2irfRiKy/n In ■ 

(16) 

with 7 the stacking fault energy, / a shape fac¬ 
tor close to unity, r c the dislocation core radius 
close to the Burgers vectors, and K an average 
value of the prefactor appearing in the elastic en¬ 
ergy. The values of these parameters for the dif¬ 
ferent vacancy loops in zirconium can be found in 
Ref. j20|. The quantity R\ is a scaling distance 
with R\ = a(-\/3/27r) 1 / 2 for the basal loops and 
R\ = y/ ac/2ir for the prismatic loops. 

As suggested in Refs. [Til . [Tl . [r0 | , the effect of hy¬ 
drogen on cluster stability can be included in a very 
simple way in the Eq. 1161 through the modification 
of the stacking fault energy, which is the leading 
term for large vacancy clusters. This is done by re¬ 
placing 7 in Eq. [THl by the modified stacking fault 
energies 7({c p }), calculated using Ea. fl5l 

In Fig. [9] the predicted formation energies of the 
loops lying in the basal and prismatic planes are dis¬ 
played in the absence of hydrogen and for a hydro¬ 
gen nominal concentration xn = 0.01 at T = 600 K. 


E 

o 

u- 



(a) Without hydrogen 



Figure 9: Formation energies of large vacancy loops as pre¬ 
dicted by continuous laws parameterized on DFT results for 
(a) pure hep Zr and (b) Zr containing hydrogen at T — 600 K 
for a nominal concentration rg = 0.01. The H vibration free 
energy is included in (b). The result obtained without this 
vibration free energy is also shown by a thin continuous line 
for the basal loops with an E fault. The insets show the sta¬ 
bility inversion at small sizes between basal and prismatic 
loops. 


The presence of hydrogen tends to increase the sta¬ 
bility of (c) loops lying in the basal planes with 
respect to (a) loops lying in the prismatic planes. 
The effect is stronger for (c) loops with an E fault, 
because of the highest amount of hydrogen captured 
by this stacking fault. These (c) loops with an E 
fault are the most stable vacancy clusters at very 
small size (loops containing less than 16 vacancies 
in pure Zr) and H segregation shifts the stability 
crossover with prismatic loops to larger sizes, an ef¬ 
fect enhanced by the contribution of H vibrations. 
The hydrogen segregation on stacking faults can 
thus explain the higher amount of (c) component 
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loops observed when Zr samples are pre-hydrided 
and then submitted to irradiation [161 Il7|] , and the 
correlated consequence on the breakaway growth at 
high irradiation dose [TiJ. 

6. Conclusion 

The influence of hydrogen in solid solution on va¬ 
cancy loops energetics has been studied in hep Zr 
using an ab initio approach, including H vibrations. 
DFT calculations performed on small vacancy clus¬ 
ters showed that hydrogen is able to attract up to 
4 vacancies. The most favorable insertion site for 
H in the vicinity of a vacancy cluster is not a T site 
but an O site with 2 or 3 vacancies as first nearest 
neighbors. Hydrogen therefore prefers to sit at the 
boundary of vacancy clusters, and not inside them. 
No strong discriminating effect of hydrogen on the 
stability of the different cluster configurations were 
found for small vacancy clusters. 

Hydrogen binding to larger vacancy clusters has 
been modeled through its interaction with the 
stacking faults. We computed, with ab initio cal¬ 
culations, the interactions between H and the pris¬ 
matic, basal E and basal Ii faults. These interac¬ 
tions are mainly attractive. Hydrogen binding is 
stronger for the basal faults, especially in the case 
of the E fault. Using these binding energies to¬ 
gether with a thermodynamic modeling of H segre¬ 
gation on the faults, we could predict the evolution 
of the stacking fault energies with the H nominal 
concentration. Including these variations into con¬ 
tinuous laws describing the energetics of vacancy 
loops, we could finally demonstrate that H stabi¬ 
lizes the vacancy loops, and that the stronger vari¬ 
ation is observed for the basal loops, i.e. the defects 
responsible for the breakaway growth observed at 
high irradiation dose. This binding of H to faulted 
vacancy loops agrees with the trapping of hydro¬ 
gen by vacancy clusters observed experimentally in 
irradiated zirconium alloys. 

Appendix A. Convergence with supercell 
size 

In Fig. I A. 101 the binding energies of H with an 
octahedral V6 cluster are displayed versus the num¬ 
ber of atoms in the supercell. The three more fa¬ 
vorable insertion sites are tested for the H atom: T 
site with one nearest neighbor HV interaction, and 
O site with 2 or 3 nearest neighbor HV interactions 
(p = 2 and 3, respectively). 



Figure A. 10: Binding energy of H with an octahedral va¬ 
cancy cluster V6 as a function of the supercell size. Three 
different insertion sites are shown for H: T site with one near¬ 
est neighbor HV interaction and O site with p = 2 or p = 3 
nearest neighbor HV interactions. 


The resulting binding energies only slightly vary 
between 100 and 200 atoms, assessing thus our 
choice of a 200 atom supercell (5 x 5 x 4) to study 
the influence of H on small vacancy clusters. 


Appendix B. Validity of the harmonic ap¬ 
proximation 


As pointed out by Christensen et al. [22], T in¬ 


sertion sites in the hep lattice form pairs of close 
sites and the migration barrier for H jumping from 
one T site to its neighbour site is small. Using the 
simulation cell containing 200 lattice sites we com¬ 
pute the H migration energies between its different 
insertion sites (Tab. El and obtain the lowest en¬ 
ergy for the T-T migration. This migration energy 
is in the range of the zero point energy of H vi¬ 
bration. One can therefore wonder if the harmonic 
approximation used to calculate the contribution of 
H vibrations to the free energy is a good approxi¬ 
mation. We examine this point more closely in this 
appendix, first for a single H atom in a perfect crys¬ 
tal, then for a H atom close to a vacancy. 


Appendix B.l. Single H atom 

We show in Fig. IB.Ill the energy profile of an 
H atom between two neighbour insertion T sites. 
Only the equilibrium configurations, corresponding 
to xh/A = 0 and 1 with A = c/4, have been relaxed. 
They are then used to obtain the other intermedi¬ 
ate configurations by linear interpolation. The ob¬ 
tained energy profile is well fitted by an order 4 
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Table B.4: H migration energies _E ml s between neighbour T 
and O sites. 



^mig 

(ev) 


unrelaxed 

relaxed 

T -o- T 

0.21 

0.13 

T —s> 0 

0.46 

0.42 

0 0 0 

0.54 

0.50 



Figure B.ll: Energy variation of a H atom migrating be¬ 
tween two T sites in a hep Zr lattice. The double well profile 
given by the solid line has been fitted to the ab initio energies 
indicated by filled circles and the corresponding quantized 
energy levels are given by horizontal solid lines. The results 
obtained for a single well in the harmonic approximation are 
shown with dashed lines. 

polynomial function. We then compute the eigen¬ 
states of this double well energy profile using the 
numerical approach described in [42[. The double 
well shape has two main effects on these eigenstates 
compared to the results obtained with the harmonic 
approximation for two single wells (Fig. IB. Ill) : it 
lowers the energy of the eigenstates, by 6meV for 
the fundamental level, and it lifts the degeneracy of 
the higher levels. 

We then compute the H vibration free energy, 
following the approach proposed by Christensen et 
al. [22j . We factorize the partition function into vi¬ 
bration modes corresponding to the H migration 
between two neighbouring T sites and vibration 
modes orthogonal to this migration direction. The 
first contribution is deduced from the eigenstates of 
the double-well energy profile ('Fig. IB.lTI) . whereas 
the harmonic approximation (Eq. [T|) is used for the 
two last vibration modes. The obtained H solution 
free energy is shown in Fig. IB. 121 where we compare 



Figure B.12: Solution free energy F (Eq. [3j of an H single 
atom in a Zr hep matrix as a function of the temperature 
calculated for its T and O insertion sites. The free energy is 
computed either within the harmonic approximation or using 
the eigenstates of the double well energy profile (Fig. I [ 1.1 I l i. 

it to the solution free energy obtained within har¬ 
monic approximation. The difference between the 
two plots is mainly a constant energy shift corre¬ 
sponding to the shift of the fundamental eigenstate 
noticed before. The linear decrease with the tem¬ 
perature of the solution free energy close to 0 K is 
well captured by the harmonic approximation, once 
the —kT In (2) has been included in the free energy 
to account for the configurational entropy associ¬ 
ated with the presence of two energy wells. The 
harmonic approximation looks therefore reasonable 
to account for H vibration in an hep Zr matrix, with 
a precision of the order of 10 meV. This precision 
being larger than the energy difference between the 
T and the O configurations (Fig. IB. 121) . a more 
precise treatment of vibrations will nevertheless be 
needed to be able to conclude on the preferential 
insertion site of H in Zr. 

Appendix B.2. H close to a vacancy 

We now examine the case where the H atom lies 
in a T site close to a vacancy. We examine two 
different configurations, both corresponding to the 
H atom and the vacancy being first nearest neigh¬ 
bours. The HV pair can either form a configuration 
b (see Fig. |T] for a sketch of the configuration), with 
both wells of the T positions corresponding to this 
b configuration, or the HV complex can form a con¬ 
figuration a, with the other well corresponding to 
configuration c. The corresponding energy profiles 
are shown in Fig. IB.131 The presence of the va- 
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Figure B.14: Binding free energy v between a vacancy 
and an H atom lying in a double-well composed of two T 
sites with a configuration either b-b or a-c. The free energy 
is computed either within the harmonic approximation or 
using the eigenstates of the double well energy profiles (Figs. 
IB. Ill and IB. 131 . 

file along the migration direction. For the asym¬ 
metric double-well energy profile corresponding to 
the transition between the configurations a and c 
of the HV complex (Fig. IB. 14b h one cannot sepa¬ 
rate the vibration mode in the migration direction 
from the modes orthogonal to this direction. An 
approximate treatment is thus needed. For such a 
double-well, the free energy is defined by 


Figure B.13: Energy variation of a H atom migrating be¬ 
tween two T sites in the vicinity of a vacancy in a hep Zr 
lattice. In a), both equilibrium positions correspond to b 
configurations (Fig. HI . In b), :r h = 0 and 1 correspond 
respectively to the a and c configurations. 


cancy leads to a decrease of the migration barrier 
between the two H equilibrium positions compared 
to the migration barrier in the bulk crystal (Fig. 
IB.Ill) . As a consequence, the effect of the potential 
anharmonicity are more pronounced in presence of 
a vacancy: the lowering of the eigenstates is more 
important and the degeneracy is lifted even for the 
fundamental eigenstate. 

We then compute the binding energy between 
H and a vacancy for these two configurations, in¬ 
cluding in Eq. [I] the vibration contribution in the 
free energy of the ZrHV and the ZrH configura¬ 
tions. Like in the previous case, we consider either 
a fully harmonic approximation, or we include the 
eigenstates corresponding to the exact energy pro- 


F(T) = —kT In 



FZ ih \ 

kT J 


- exp 


A E + F' r 

W 


vib \ 


(B.l) 


where the vibration contributions F!^ lb and F^ lh are 
given by Eq. [T] for both configurations a and c, and 
A E = E c — E a is the energy difference between 
both configurations. The free energy deduced from 
the exact eigenstates of the ID double-well (Fig. 
IB. ldll can be added to this free energy, after having 
withdrawn the corresponding contribution in the 
harmonic approximation: 


F${T) = -kT In 


1 

2 sinh (Hui^/kT) 
exp (—A E/kT) 
^2 sinh (huj^/kT) 


(B.2) 


with wf and the pulsations of the vibration 
modes in the migration direction. 
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The binding free energies are shown in Fig. IB . 141 
We recover a difference of the order of lOmeV be¬ 
tween the result of the harmonic approximation and 
the more exact treatment of H vibrations. This dif¬ 
ference is much smaller than the absolute values of 
the binding energies and does not modify the rela¬ 
tive stability between different configurations. The 
harmonic approximation appears clearly to be well 
suited to study H interaction with vacancies and a 
more precise treatment of H vibrations is thus un¬ 
necessary. 
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